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Possible universal dynamics of a many-body system far from thermal equilibrium are explored. A focus is set 
on meta-stable non-thermal states exhibiting critical properties such as self-similarity and independence of 
the details of how the respective state has been reached. It is proposed that universal dynamics far from 
equilibrium can be tuned to exhibit a dynamical transition where these critical properties change 
qualitatively. This is demonstrated for the case of a superfluid two-component Bose gas exhibiting different 
types of long-lived but non-thermal critical order. Scaling exponents controlled by the ratio of 
experimentally tuneable coupling parameters offer themselves as natural smoking guns. The results shed 
light on the wealth of universal phenomena expected to exist in the far-from-equilibrium realm. 

The concept of universality has been extremely successful in classifying and characterising equilibrium states 
of matter. For example, there are different types of order in a magnetic material separated by a second- order 
phase transition at which the relevant physical properties become independent of the microscopic details of 
the system. This constitutes universality and allows to characterise an extensive range of different phenomena in 
terms of just a few classes governed by the same critical properties. In view of the intensifying discussion on the 
dynamics of many-body systems it is a pressing question whether also far away from the thermal limit the 
character of dynamical evolution can become independent of the microscopic details. For a closed system this 
would imply that in approaching a critical configuration the evolution must become independent of the particular 
initial state the system has started from and critical slowing down in the actual time evolution is observed. As a 
consequence, different types of dynamical evolution could be distinguished by means of universality classes. 

In this article, we demonstrate that such universal dynamics far from thermal equilibrium states is indeed 
possible. We show that there is a dynamical transition between different types of universal dynamics. After having 
quenched a two-component ultracold Bose gas we follow the ensuing evolution of the closed system and observe 
transient universal behaviour. In particular, we identify a tunable external parameter that determines the type of 
transient non- equilibrium order. In our setup non- equilibrium order is constituted by the appearance of spatial 
patterns including quasi-topological defects like vortices and skyrmions which are created through instabilities 
and which decay only on very large time scales 1 . We find the dynamics of diluting defects to constitute a separate 
dynamical critical phenomenon far from thermal equilibrium and therefore clearly distinct from the known 
equilibrium critical points 2 ' 3 . Large-scale correlations in the system are universal in the sense that they are fixed by 
the type of defect, but are completely insensitive to the specific positions and velocities of the defects. 
Distinguishable types of defects are produced for different values of the tuning parameter. Hence, the dynamics 
can be tuned to a transition between different metastable non -equilibrium ordered states, see Figure 1. Extending 
upon the presented results will allow to learn about the universal properties of classes of models in the vicinity of 
attractors, non-thermal fixed points and other critical submanifolds of the greater space of all possible non- 
equilibrium configurations 4 " 7 . Knowing such universal properties, predictions for the behavior of very different 
physical systems can be obtained on the basis of comparatively few exemplary measurements. When looking at 
fundamental science applications, universality classes of far-from-equilibrium criticality can link the dynamics, 
which we propose to be observable in ultracold atomic or exciton-polariton gases, with phenomena as different as, 
e.g., magnon gases in solids, quark-gluon-plasma dynamics in heavy-ion collisions, or reheating after early- 
universe inflation. Technological applications would not stay away. 

The two-component Bose systems we consider are described by the Hamiltonian density (h = 1) 



(i) 



Here m is the mass of the atoms and a = g l2 /g the ratio between the inter-species coupling g 12 and the intra- 
species interaction constant g. The latter as well as the mass are chosen to be the same for both components and 
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Thermal states 



Figure 1 | Illustration of the dynamical evolution of a two-component Bose gas, with a focus on the possible transition between far-from-equilibrium 
transient states. The set in the bottom left corner represents an ensemble of out- of- equilibrium initial states. Amongst these, a subset subjected to a 
dynamical instability is marked in white. It is chosen such that for all values of the the ratio a of inter- to intra-species couplings the system belongs to the 
unstable subset. Subsequently, the time evolution of the closed system leads to non- and quasi-topological defect formation. The type of defects created 
determines the scaling in the low-momentum regime of the particle spectrum of the gas. Our results reveal three different scalings corresponding to three 
types of defect configurations, and a dynamical transition between the far-from-equilibrium universal states. 



the sum over the field index j e {1, 2} is implied. The considered 
systems allow for good experimental control and have been studied 
intensely 8 " 11 . In experiment, a can be varied by means of a Feshbach 
resonance. Two different ground states exist, depending on the value 
of the parameter a 1213 . In the immiscible regime, a > 1, the inter- 
species interaction energy is greater than the intra-species interac- 
tions. Hence, in the ground state, the spatial overlap of the 
components is minimised through domain formation and spatial 
separation of the two components. On the contrary, for a < 1, the 
two components become miscible and, in the ground state, uni- 
formly distributed over the whole volume. We make use of a to 
change the properties of the system in the yet unexplored region of 
non-equilibrium quasi- stationary states. A focus is set on long-lived 
states with non-/quasi-topological defects including domain walls 812 , 
single- species vortices, and skyrmions 14 in the coupled system. 

Results 

Since the dynamics we are interested in exclusively affects the low- 
momentum, strongly populated field modes, we employ the so-called 
classical field method which yields, within numerical accuracy, exact 
results for the time- evolving observables. In order to discuss the 
contribution of the different (quasi-)topological configurations in 
detail, we make use of the spin representation of the two -species 
fluid. In this representation the S z - component of the spin vector 
simply refers to the density imbalance between the two species. 
This is an essential observable for the detection of domain walls. 
Furthermore, we make use of a hydrodynamic decomposition of 
the superfluid flow which allows for the detection of vortex contribu- 
tions to the spectra. For details on numerical parameters and decom- 
position methods we refer the reader to the Methods section. In 
Figure 2, we unravel the time evolution of a two-component Bose 
gas as outlined in Figure 1. The parameters of our simulation are 
chosen such that the final states are close to the ground state of the 
system. To clarify the type of non-equilibrium order during the 
intermediate stages of the evolution, we show the imbalance between 
the two components S z as well as the momentum distribution n(k) for 
three characteristic times. The initial time, t h marks the stage of an 
isotropic momentum distribution which is overpopulated within an 
intermediate range of momenta, as compared to the ensuing equi- 
librium distribution. It is marked by a strong fall-off at large 
momenta. In our driving scheme this state is reached in the wake 



of an instability. For the immiscible case, a > 1, we use a modula- 
tional instability in which small low-momentum fluctuations in the 
polarisation are amplified to macroscopic spin domains. In the mis- 
cible regime, a < 1, we invoke a counter- superflow instability 1516 by 
choosing oppositely directed flow vectors for the two field compo- 
nents. Above a critical relative velocity, momentum exchange 
between the two fluids causes the superflow to decay and spin 
domains to form. Subsequently, overpopulation at intermediate 
momenta is encountered for all values of the parameter a. Its micro- 
scopic origin is seen in the coloured distributions of the density 
imbalance S z in Figure 2 which show strong fluctuations and the 
onset of domain formation. We estimate t\ ~ 10/<x>i, with ojj being 
the energy of the fastest growing unstable mode. Note, that the 
dynamics up to t Y is different for different a and involves for example 
isotropisation in case of a < 0. For details see section A in the 
Supplementary Material. The further dynamical evolution of over- 
populated states involves particle transport towards small momenta 
and energy transport towards high momenta 1718 . In Figure 2 this is 
best observed by comparing the high momentum tails of the 
momentum distributions at t\ and £ N tfp- Particle transport towards 
low momenta eventually fills the zero -mode (not seen), a process 
most important in the miscible regime a = 0.8. The key insight 
gained from the spectra concerns the development of infrared (IR) 
power laws in the momentum distribution n(k) ~ k~ c at 
f NTFP ~ 10 2 /cl>i. These power laws depend on the external parameter 
a: C ~ 3.5 for a > 1, £ ~ 3 for a = 1 and £ ~ 4 for a < 1. Also, the 
pattern of imbalance fluctuations in the system depends on the inter- 
action strength. 

We highlight the most striking observations at late times before we 
carry on developing a detailed understanding of the relation between 
the dominant fluctuations and the observed power laws. In the 
immiscible regime the system consists of few large domains with 
additional inclusions of small point-like domains. As compared to 
the initial time domains have considerably grown. A special situation 
is encountered at the transition point (a = 1), where domain-like 
structures persist to extremely long times. This is remarkable since in 
this regime domains are not energetically favourable as compared to 
overlapping particle densities. For a = 0.8, imbalance fluctuations 
have decayed up to few small areas of strong imbalance. Let us finally 
look at the largest computed time £ F — 10 3 /<x>i. For a > 1, two 
domains of equal size remain, which reflects the immiscibility in 



SCIENTIFIC REPORTS | 3 : 2394 | DOI: 1 0.1 038/srep02394 



2 




t\ ^ntfp time 



Figure 2 | Time evolution of the system in three different parameter regimes. Three snapshots (columns) of the evolving two-component Bose gas for 
three different values of the coupling ratio a (rows) are shown. In each case we show the spatial distribution of density imbalance S z as well as the 
momentum distribution n(k) as a function of the momentum modulus k of particles in the gas. Tick labels for the spectra graphs are the same as in 
Figure 3. Until ti the dynamics is characterised by the instability and isotropisation of density fluctuations. At *ntfp = 10 *i> one observes the development 
of different power laws in the momentum distributions. They remain metastable for a long time beyond t F = 100 t\ and reveal different non-thermal fixed 
points. The change of power laws at a = 1 indicates a dynamical transition. 



the ground state. The small-scale domain-fluctuations have consid- 
erably reduced in number. At the transition point, we observe the 
persistence of domain-like structures which we attribute to a diver- 
ging time scale t of domain decay as a — » 1 from below. For details see 
section B in the Supplementary Material. For a < 1, small long-lived 
imbalance fluctuations remain, whereas the background tends to 
become very smooth, S z ~ 0. Finally, we investigate the microscopic 
origin of the scaling found in the momentum distributions. The 
bimodal power laws in these distributions are signatures of the sys- 
tem having approached a non-thermal fixed point 4 " 7 . At long times 
they become more pronounced in all three cases, signalling critical 
slowing down of the dynamic evolution. In Figure 3, the result of a 
decomposition of the momentum distribution according to spin and 
fluid degrees of freedom at t = £ntfp is presented, see the Methods 
section for details on the decomposition. We show spin- spin correla- 
tion functions as well as correlations in the incompressible velocity as 
a function of momentum. The crucial point is that they explain the 
scaling properties of the total momentum distribution and allow for 
an interpretation in terms of specific defect configurations, as shown 
in Figure 2. This completes our understanding of transient non- 
equilibrium order in a two-component ultracold Bose gas. 

We remark that, for all a, compressible and quantum pressure 
excitations, which are not shown, dominate the spectrum in the 



ultraviolet (UV) regime and follow a thermal Rayleigh-Jeans distri- 
bution n ~ /c" 2 , while they give a negligible contribution towards 
lower momenta. 

In the immiscible regime, see Figure 3 (left), the main contribution 
to the spectrum in the IR region is provided by the incompressible 
component n b corresponding to flow orientations transverse to the 
direction in which the flow velocity varies. Although we are dealing 
with a multi- component gas this feature is similar to superfluid 
turbulent flow in a single- component Bose gas. Thereby the incom- 
pressible spectrum shows n { ~ k~ 4 scaling over approximately one 
decade which is generated by coherent vortical flows w,- around 
topological defects 1718 . The cores of these vortex-like structures can 
be seen in the spin imbalance, see Figure 2, since they are filled with 
particles of the other species and thus are of the skyrmion type. They 
are created during the merging process of domains and, persisting 
due to their topological nature, give the main contribution to the 
incompressible component spectrum. However, spin excitations n z s 
overtake in an intermediate momentum region, showing a scaling of 
n z s ~k~ 3 . Looking more closely, one observes that the scaling beha- 
viour terminates in the IR at a scale 7i/L D given by the mean domain 
size L D , while the cut-off in the UV at n/£ s is set by the width of the 
domain walls, i.e., the spin healing length £ 5 = f (l/|l — a|) 1/2 . Since 
the two contributions n z s and n t are of comparable magnitude within 




Figure 3 | Decomposition of correlation functions of the two-component Bose gas. The three different graphs refer to the three different non- 
equilibrium ordered states parametrised by the coupling ratio a. Shown are the z-spin component n z s > the sum of x- and y- spin components n* y and the 
incompressible component « f at time Jntfp- F° r details of the decomposition procedure see the Methods section. 
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an intermediate momentum range the sum of all contributions, giv- 
ing the full spectrum n(k), appears to follow the power law n ~ /c" 3 5 
in the IR, as discussed above. 

The situation on the other side of the transition can also be cla- 
rified. When a falls below 1 domains are energetically suppressed. 
Thus, vortices dominate the non-thermal fixed point and induce the 
characteristic scaling n ~ n { ~ k~ 4 , see Figure 3 (right). Here, it is also 
visible in the scaling of relative phase fluctuations n* y , which is 
related to constant particle densities on scales larger than £ as com- 
pared to the skyrmionic case above. 

The picture changes dramatically on the transition point, a = 1. As 
shown in Figure 3 (middle), vortical flow is much less important in 
this regime. The momentum distribution is dominated by spin fluc- 
tuations scaling as n ~ n z s ~ n s y ~ k~ 3 . We attribute this feature to the 
increased dimensionality of the vacuum manifold for a = 1 which 
removes the topological protection of vortices. The stability of 
domain walls at the transition point takes over but is non-topological 
in nature. Instead, conservation laws restrict the decay of this par- 
ticular type of defect 19 . This argument can be made even more trans- 
parent by studying the critical dynamics at the transition from a < 1 
towards a = 1, where two non-equilibrium ordered states meet each 
other, vortices and domains. Details on the connection between 
defect structures and scaling exponents of the momentum distri- 
bution can be found in the Methods section. 

Discussion 

We conclude that the appearance of (non-) topological defects deter- 
mining bulk features of correlation functions in far-from-equilib- 
rium situations is very general and strongly relevant for 
intermediate-time universal dynamics. It is easily imaginable that 
multi- component field theories with more than two components 
show similar behavior to the one outlined here. We point out the 
possibility of a universal duality between diluting defects and inverse 
particle cascades. This requires the generation of (quasi- topological 
configurations far from thermal equilibrium and their slow decay, 
going together with an increase of coherence and defect separation 20 . 
Under these conditions, an inverse particle cascade is generated, and 
the associated power-laws can be found from the scaling properties 
of the respective single defect. 

The cascades exhibit a characteristic non-thermal infrared power- 
law behaviour of correlation functions, e.g., as demonstrated in the 
paper, the momentum spectrum n(k) ~ k~ c , £ ~ 3 . . . 4. This scaling 
is distinctly different from the situation of a thermal gas and reflects a 
different relation between spectral and statistical correlations far 
from equilibrium. The scaling behaviour which is evolving further 
towards the infrared while the defects are being diluted marks the 
presence of a non-thermal fixed point. In the vicinity of this point, 
the time evolution becomes universal due to self- similarity in the 
scaling regime. We point out that universality, here, is not only 
related to the scaling in the spatial coordinates but is also present 
in the time dimension, manifesting itself e.g. in the cascading char- 
acter of the scaling regime which is progressively building up coher- 
ent population in the infrared due to the dilution process of the 
topological defects. Furthermore, the time evolution of correlation 
functions can show scaling behaviour. We emphasise that adiabatic 
near- equilibrium descriptions of dynamical critical phenomena do 
not apply to the strong quench scenarios studied here, such that in 
principle new scaling exponents should appear which mark new 
universality classes of far-from-equilibrium critical dynamics, the 
determination of which is beyond the scope of this work. 

We furthermore conclude that it is possible to establish the notion 
of a dynamical transition between different far-from-equilibrium 
ordered states in a two-component Bose gas. Tuning a system para- 
meter (a) through a point where the symmetry properties of the 
system change allows the universal far-from-equilibrium dynamics, 
e.g., its scaling properties, change qualitatively. Hence, in what we call 



a dynamical transition, the way how non-thermal fixed points are 
approached during the equilibration evolution changes qualitatively. 
Determining how such a dynamical transition is related to a corres- 
ponding equilibrium transition, being beyond the scope of this work, 
poses an interesting question for the development of a theory of far- 
from-equilibrium critical phenomena. 

A variety of (quasi-)topological excitations are known to exist in 
multi- component fields 119 , examples are monopoles in gauge fields 21 
and exotic magnets 22 , as well as skyrmions in Bose-Einstein conden- 
sates 14,23 and liquid crystals 24 . New interesting features that are read- 
ily accessible in experiment are expected for ultracold spinor 
gases 25 " 29 . The transition between different types of transient non- 
equilibrium order can be controlled by changing the symmetry prop- 
erties of the Hamiltonian and thus topology and local conservation 
laws of the system. This offers interesting prospects for far-from- 
equilibrium dynamical transitions in very different areas of physics. 

Methods 

Classical field equations. Since the dynamics we are interested in exclusively affects 
the low-momentum, strongly populated field modes, we employ the so-called 
classical field method which yields, within numerical accuracy, exact results for the 
time-evolving observables 30,31 . For this, initial field configurations $i y2 {x, t 0 ) are 
sampled from Gaussian probability distributions and then propagated according to 
the classical equations of motion. At the end of the time evolution correlation 
functions are obtained from ensemble averages over the set of sampled paths. The 
classical equations of motion derived from the Hamiltonian density (1) of the 
interacting two-component Bose gas are 

iW, = -ivVi+*(|*,| 2 + a#2|%, (2a) 

i^2=-^V 2 0 2 +g(|^ 2 | 2 + ^ 1 | 2 )^ 2 . (2b) 

The momentum distribution which is shown in various graphs is defined as 

n(fc,0=|dQ fc (^(k,%.(k,0), (3) 

where J dQk denotes the angular average in two-dimensional momentum space. 

We have applied the following rescalings to obtain dimensionless quantities in the 
system of coupled Gross-Pitaevskii equations (2), using the lattice constant a s of the 
computational grid: (j>a s — > 0, mg^> gand tj (ma]) -*t. Here n = (N 1 + N 2 )/L 2 is the 
mean total particle density onaJV s X N s grid of linear size L = N s a s , with N s = 1024. 
Conversion to physical time scales appropriate for a specific experiment, with, e.g., 
87 Rb atoms, can be achieved by noting that the healing length £ = (2mgn)~ m is given 
by £ = 4a s for our system. Giving an estimate, e.g., for a physical healing length 
£~ 1 fim, our lattice unit corresponds to a s ~0.25 fim which yields a time unit of 
m^b^/h— 10 ~ 4 s. Hence, the characteristic time scale is 1 / coi ~ 10 ~ 3 ... 10 ~ 2 s 
which yields a realistic time scale for the appearance of the non-thermal fixed point at 
£ N tfp ^ 0.1 ... 1 s. 

Spin-fluid representation. Although we always simulate the full set of equations (2) 
it is instructive to study the degrees of freedom which describe the relative evolution 
of the two components. Writing the fields in the polar representation 
h = \/ r Pi ex V (ify) these are given by the local phase difference 6 r = 6 1 — 6 2 and by the 
local density difference pi — p 2 . The Pauli matrices o a serve to define the Sch winger 
representation S a = ^\o a i ^ (sum over repeated indices implied) of angular 
momentum operators. This results in a three-component vector of (pseudo-)spin 
densities S a , a e {x, y, z} 

S x = 2 A /p 1 p 2 cos(0 r ), (4a) 
Sy=-2y/p^sm{6 r ), (4b) 

S z = P l ~P 2 , (4c) 

where the modulus corresponds to the total density |S| = p x + p 2 = p T . For 
convenience, we apply the redefinition S a — > pjS a such that |S| = 1. Then, the 
Hamiltonian density H of Eq. (1) can be rewritten as 



which shows the influence of the relative degrees of freedom and their coupling to the 
global ones 32 . The quantity j T = pi V^i + p 2 ^(p 2 is the conserved total particle current 
associated with a global (7(1) symmetry of Eq. (1) associated with the global shift of 
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the total phase 6 T = 9 X + 6 2 . Thus ) T can not be expressed using just the spin densities 
but contains also the remaining fourth degree of freedom, the total phase 6 T . We 
obtain 

i - - Pt S70 t + - PtSZ 91 (S? VS x - S x VSH . (6) 
Jr 2^ r 2[l-(S z ) 2 ] V ' V 7 

According to the representation (5) of the Hamiltonian density the two-component 
Bose gas can be described by a spinor field S which is coupled to a hydrodynamic fluid 
with density p T and (conserved) quasi particle current j T . For a fluid at rest, i.e., p T = 
const and ) T = 0, the spin system thereby assumes the form of a classical non-linear 
sigma model (NLaM) with a mass term gp\ /4( 1 — a) [(S x ) -\-(S y ) ] , whereas in 
general a current j T ^ 0 leads to a highly non-trivial coupling between internal and 
hydrodynamic degrees of freedom. 

Hydrodynamic decomposition. In order to gain further insight into the interplay 
between the domain structure and other excitations we decompose the energy 
spectrum according to Refs. 18,33. The results of this decomposition at the non- 
thermal fixed point were discussed in Figure 3 of the main text. In a density-phase 
representation, the kinetic part of the Hamiltonian density can be written in the 
following way: 

n kin = vtyvtj = |Vp7| 2 + ^ ws a ws a + |w| 2 . (7) 

The velocity field w is defined via the total particle current y/p^w = j r , similar to the 
convenient choice for the one-component case. In this decomposition the first and 
the last terms are the quantum-pressure component and the classical hydrodynamic 
component, respectively. In contrast to a single-component fluid the second term of 
Eq. (7) adds a pressure-like contribution to the kinetic energy which is produced by 
internal excitations only. In addition, the velocity field w and thus the corresponding 
part of the kinetic energy can be decomposed in a compressible and an incompressible 
part, w = w c + w,-, with V X w c = 0 and V • w,- = 0, such that wave-like and vortical 
excitations appear in different parts of the decomposition, see Refs. 33. Based on Eq. 
(7), radial energy spectra in momentum space corresponding to the discussed energy 
parts can be defined as 

E 5 (k)= 1 -^dQ k (\w 5 (k)\ 2 ), 5e{q,c,i} (8a) 

E s (k)= l -^dQ k {w a s (k)-w a s (k)). (8b) 

Here we have introduced additional velocities = Vy^Pr an d = ^/p^S/S a /2 for 
the sake of a closed representation. Finally, the energy spectra can be converted to 
occupation number spectra defined in Eq. (3) by multiplication with a factor k~ 2 , 
n s (k) = k- 2 E s (k), S e= {q, c, i, s} 18 . 

Domain and vortex spectra. In Figure 3 of the main text we discussed the 
microscopic origin of characteristic power-laws in the momentum distribution n(k) 
by decomposing the spectra in different components. This was successful since 
signatures of different defects show up in specific components of the spectra. Here we 
make this correspondence explicit for the case of vortices and domain walls. 

Under the assumption of spatially constant total particle density the spectrum n s 
can be related to the Fourier transform of the correlation function of the angle- 
averaged spin order parameter S(k) = § dQk(S a ( — k)S fl (k)), 

k 2 n s =^ [dQ fc (^[VS fl *]^[VS fl ]> 

f (9) 
= ^k 2 jdQ k (S a (-k)S a (k)), 

with T denoting the Fourier transform, and therefore n s (k) = p T <S(k) /2. Hence, in a 
momentum regime 7i/L D <C /c <C ft/ £ s where the the scaling behaviour is dominated by 
a single domain wall, n s ~S~k~ 3 follows from an ansatz in terms of the Heaviside 
function S z = ± (1 — 29{x)), S x = S? = 0. This feature is similar to the scaling induced 
by solitons in one- dimensional Bose gases, where the phase jump occurs in the 
bosonic field and induces a scaling n 1D ~ k~ 2 , see Ref. 34. 

Similarly, the appearance of quantised vortices induces scaling behaviour of the 
incompressible energy n { . On scales considerably smaller than the mean distance 
between vortices, the effective velocity field |w,-(r)| associated with the rotational flow 
decays as II r with growing distance r from the nearest vortex core. As a result, the 
angle-averaged incompressible kinetic energy distribution ~ |w 2 \/l gives rise to the 
IR momentum spectrum n^k) ~ k~ A , whereas compressible excitations vanish. 
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